Теория вероятностей
Метод Монте-Карло
Сущность метода Монте-Карло
Метод Монте-Карло основан на использовании случайных чисел. Он применяется во многих областях науки и техники (статистическая физика, физика твердого тела, физика плазмы, астрофизика, теория массового обслуживания, теория игр, теория управления и др.) Метод Монте-Карло используют для вычисления интегралов, в особенности многомерных, для решения систем уравнений, для решения задач оптимизации, для исследования различных сложных систем (экономических, биологических и т.д.). Это один из наиболее мощных и универсальных методов решения сложных задач из самых разных областей знаний.
Название метода связано с названием района Монте-Карло, знаменитого своими казино, где играют в рулетку — одно из простейших устройств для получения случайных чисел.
Сущность метода Монте-Карло состоит в следующем: требуется найти значение некоторой изучаемой величины. Для этого выбирают такую случайную величину , математическое ожидание которой равно . Производят испытаний, в результате которых получают возможных значений случайной величины , вычисляют их среднее арифметическое
и принимают в качестве оценки (приближенного значения) искомого числа :
Особенности метода Монте-Карло:
- Простая структура вычислительного алгоритма. Как правило, некоторое случайное испытание повторяется раз, и результаты всех опытов усредняются. Поэтому иногда метод Монте-Карло называют методом статистических испытаний.
- Погрешность вычислений пропорциональна , где — число испытаний. Например, для того чтобы уменьшить погрешность в раз (т.е. чтобы получить в ответе еще один верный десятичный знак), нужно увеличить (т.е. объем вычислений) в 100 раз. Для достижения приемлемой точности число испытаний должно быть достаточно большим.
Генерация случайных чисел в Engee
Компьютерные алгоритмы генерируют не истинно случайные числа, имеющие бесконечное число десятичных знаков, а псевдослучайные числа, имеющие конечное число знаков.
Генерация равномерно распределенных случайных чисел
Функция rand(dims) генерирует массив размером dims со случайными, независимо, одинаково и равномерно распределенными числами в интервале . При отсутствии аргументов функция rand() генерирует одно случайное равномерно распределенное число в интервале .
rand()
rand(10)
rand(4, 5)
Функция rand(a:h:b, dims) генерирует массив размером dims со случайными равномерно распределенными числами, принимающими возможные значения в диапазоне от a до b с шагом h. Если шаг h не задан, то он принимается равным . Если параметр dims не задан, то генерируется одно число.
rand(1:10)
rand(0:0.01:10, 10)
rand(60:0.1:70, (3, 3))
Функция rand([a,b,c,...], dims) генерирует массив размером dims со случайными равномерно распределенными числами, принимающими возможные значения только из множества элементов вектора [a,b,c,...].
Например, сгенерируем равномерно распределенных случайных чисел, которые могут принимают значения , или :
rand([1, 2, 5], 20)
Для генерации случайного числа, имеющего равномерное распределение и принимающего непрерывные возможные значения в интервале , можно использовать выражение a+(b-a)*rand().
Например, сгенерируем случайное число, равномерно распределенное в интервале :
3 + (8-3) * rand()
Генерация нормально распределенных случайных чисел
Функция randn(dims) генерирует массив размером dims со случайными, независимо и одинаково распределенными числами, имеющими стандартное нормальное распределение (с математическим ожиданием, равным нулю, и со средним квадратическим отклонением, равным ). При отсутствии аргументов функция randn() генерирует одно случайное число, имеющее стандартное нормальное распределение.
randn()
randn(10)
randn(6, 4)
Для генерации случайного числа, распределенного нормально с произвольными значениями математического ожидания a и среднего квадратического отклонения sigma, можно использовать выражение sigma * randn() + a.
Например, сгенерируем случайное число, распределенное нормально с параметрами и :
3 * randn() + 5
✏️Задание 1
С помощью функций rand и randn сгенерируйте:
- матрицу размером случайных чисел, равномерно распределенных в интервале ;
- одномерный массив из равномерно распределенных случайных чисел, принимающих значения в диапазоне от до с шагом ;
- случайное число, равномерно распределенное в интервале ;
- одномерный массив из случайных чисел, имеющих стандартное нормальное распределение;
- случайное число, распределенное нормально с параметрами и .
Решение
rand(5,5)
rand(0:10:150, 12)
100 + (200-100) * rand()
randn(15)
2 * randn() + 10
Генерация случайных чисел с произвольным законом распределения
Функция rand(p, dims) генерирует массив размером dims со случайными числами, распределенными по заданному закону p. В предыдущих разделах мы рассматривали примеры такой генерации для биномиального, равномерного, показательного, нормального распределения и для распределения хи-квадрат.
В документации к библиотеке Distributions.jl вы можете найти больше функций для генерации различных законов распределения.
Рассмотрим теперь несколько примеров применения метода Монте-Карло для решения математических и физических задач. В этих примерах мы будем использовать только равномерно распределенные случайные числа.
Вычисление значения числа
Построим в прямоугольной системе координат квадрат со стороной . Стороны квадрата расположим параллельно осям координат, а левую нижнюю вершину поместим в начало координат. Затем построим четверть окружности с центром в начале координат и радиусом .
Площадь квадрата: , площадь четверти круга: . Выразим число через отношение этих площадей:
Будем случайным образом помещать точки внутрь квадрата. Координаты и этих точек являются равномерно распределенными случайными числами в интервале . Отношение числа точек, попавших внутрь четверти круга, к общему числу поставленных точек, приблизительно равно отношению площадей соответствующих фигур:
Подставляя это соотношение в формулу для вычисления числа , получим:
Чем больше мы поставим точек, тем точнее будет это равенство.
Реализуем описанный алгоритм в следующем скрипте, поставив миллион точек:
N1 = 1e6 # общее число точек
N2 = 0 # начальное значение числа точек, попавших внутрь четверти круга
for i = 1:N1 # повторяем испытания N1 раз
x = rand() # выбираем случайным образом координату x точки
y = rand() # выбираем случайным образом координату y точки
if x^2 + y^2 <= 1 # проверяем, попала ли точка внутрь четверти круга
N2 = N2 + 1 # увеличиваем на 1 число точек, попавших внутрь четверти круга
end
end
pi_approx = 4 * N2 / N1 # вычисляем приближенное значение числа пи
Найдем абсолютную погрешность полученного результата, вычислив модуль разности между приближенным и точным значениями числа :
abs(pi_approx - pi)
Вычисление определенных интегралов
Приближенное выражение для вычисления определенного интеграла от функции в пределах от до :
где — случайные значения аргумента , равномерно распределенные на отрезке . Чем больше , тем точнее будет это равенство.
Пример. Вычислим методом Монте-Карло определенный интеграл
взяв миллион точек:
a = 0 # нижний предел интегрирования
b = pi # верхний предел интегрирования
N = 1e6 # число точек
S = 0 # начальное значение суммы
for i = 1:N # повторяем испытания N раз
x = a + (b-a) * rand() # выбираем случайную точку на отрезке [a, b]
S = S + sin(x) # вычисляем значение суммы
end
I = (b-a) / N * S # вычисляем приближенное значение интеграла
Найдем абсолютную погрешность полученного результата, вычислив модуль разности между приближенным и точным значениями интеграла. Точное значение этого интеграла равно:
abs(I-2)
✏️Задание 2
Вычислите методом Монте-Карло определенный интеграл
взяв миллион точек.
Напомним, что значение функции вычисляется в Engee как exp(x).
Решение
a = 1
b = 2
N = 1e6
S = 0
for i = 1:N
x = a + (b-a) * rand()
S = S + exp(-x^2)
end
I = (b-a) / N * S
Аналогичным образом метод Монте-Карло применяется для вычисления многомерных (двойных, тройных и т.д.) интегралов. В этом случае он имеет значительные преимущества перед другими численными методами интегрирования.
Случайные блуждания
Одномерное случайное блуждание
Случайное блуждание — случайный процесс, который описывает путь, состоящий из последовательности случайных шагов в каком-нибудь пространстве.
Модели случайного блуждания получили широкое распространение в информатике, физике, химии, биологии, экономике и социологии. Классические приложения этой модели: диффузия молекулы в газе, броуновское движение частиц в жидкости, моделирование полимерных цепочек.
Начнем с простейшей модели — одномерного случайного блуждания. Пусть частица начинает двигаться вдоль оси Ox из точки . Все шаги имеют одинаковую длину, равную . На каждом шаге по времени частица делает шаг вправо с вероятностью или шаг влево с вероятностью . Направление каждого шага не зависит от направлений предыдущих шагов.
Промоделируем одномерное случайное блуждание в течение шагов по времени и построим график зависимости координаты частицы от времени . Примем вероятности шага вправо и шага влево равными .
using Plots
N = 10000 # число шагов по времени
t = 1:N # вектор отсчетов времени
x = Vector(1:N)
x[1] = 0 # начальная координата x
for i = 2:N # повторяем испытания N раз
x[i] = x[i-1] + rand([-1, 1]) # на i-м шаге выбираем новое значение координаты x (к предыдущему значению прибавляем либо -1, либо +1)
end
plot(t, x) # строим график x(t)
Представляет интерес средний квадрат смещения частицы за шагов , который пропорционален физической величине, называемой коэффициентом самодиффузии. Рассчитаем среднюю величину квадрата смещения, усреднив величину по реализаций случайного блуждания по шагов каждое:
N = 100 # число шагов по времени
M = 100000 # число реализаций случайного блуждания
x = Vector(1:N)
function random_walk() # одна реализация случайного блуждания
x[1] = 0
for i = 2:N
x[i] = x[i-1] + rand([-1, 1])
end
return x[N]^2 # вычисляем квадрат смещения на N-м шаге по времени
end
S = 0 # усредняем по M реализациям
for k = 1:M
S = S + random_walk()
end
D = S/M # выводим результат на экран
Несмотря на стохастический характер блужданий, при разных прогонах этого скрипта величина колеблется в малых пределах около своего среднего значения.
Двумерное случайное блуждание на решетке
Рассмотрим теперь случайное блуждание на двумерной решетке. В случае квадратной решетки движение состоит из шагов одинаковой длины в четырех направлениях в один из соседних узлов решетки (влево, вправо, вниз, вверх). Для простоты примем вероятности выбора каждого направления одинаковыми.
Промоделируем случайное блуждание на квадратной решетке в течение шагов по времени и построим график зависимости координаты частицы от координаты (траекторию частицы). Движение начнем из точки с координатами .
using Plots
N = 10000 # число шагов
x = Vector(1:N)
y = Vector(1:N)
x[1] = y[1] = 0 # начальные координаты x и y
for i = 2:N # повторяем испытания N раз
z = rand(1:4) # выбираем одно из 4 направлений движения
if z == 1 # шаг влево
x[i] = x[i-1] - 1
y[i] = y[i-1]
elseif z == 2 # шаг вправо
x[i] = x[i-1] + 1
y[i] = y[i-1]
elseif z == 3 # шаг вниз
x[i] = x[i-1]
y[i] = y[i-1] - 1
else # шаг вверх
x[i] = x[i-1]
y[i] = y[i-1] + 1
end
end
plot(x, y) # строим график y(x)
Случайное падение дождевой капли
Рассмотрим частный случай двумерного случайного блуждания, при котором вероятность шага вниз больше вероятности шага вверх. Для простоты примем, что на каждом шаге координата частицы уменьшается на единицу с вероятностью, равной , а координата изменяется случайным образом, как в модели одномерного случайного блуждания.
Такое случайное блуждание можно рассматривать в качестве идеализированной модели падения дождевой капли при наличии случайных порывов ветра.
✏️Задание 3
Реализуйте в Engee алгоритм случайного падения дождевой капли в течение шагов. Координаты начальной точки возьмите равными . Особенности алгоритма:
- Задайте начальное значение координаты равным , а начальное значение координаты равным ;
- Выражение для вычисления новых значений координаты задайте таким же, как в случае одномерного случайного блуждания;
- Значение координаты на каждом шаге должно уменьшаться на единицу.
Решение
using Plots
N = 10000
x = Vector(1:N)
y = Vector(1:N)
x[1] = 0
y[1] = N
for i = 2:N
x[i] = x[i-1] + rand([-1, 1])
y[i] = y[i-1] - 1
end
plot(x, y)